#!/bin/bash
#
# Script to visualize the final temperature profile. If your system size differs from
# the current, one should revise the script accordingly.
#
# Extract box length along x direction from log file.
lx=`grep -B1 'Loop time of' log.lammps |head -1|awk '{print $8}'`
ke2T=7735.925 # 2*eV/(3*kB), to convert ke to K

# Extract the final temperature profile 
tail -200 temp.profile > final.profile

# Prepare for the gnuplot script to visualize the temperature profile.
cat > .plot.script << EOF
set term png
set out 'temp_profile.png'

set xlabel 'x (Angstrom)'
set ylabel 'Temperature (K)'

set xr [0:${lx}]
set yr [500:2000]

unset key

f(x) = a * x + b
g(x) = c * x + d

fit [ 9:35] f(x) 'final.profile' u (\$2*${lx}):(\$4*${ke2T}) via a, b
fit [48:71] g(x) 'final.profile' u (\$2*${lx}):(\$4*${ke2T}) via c, d
T = a * (b-d)/(c-a) + b

lbl_text = sprintf("T_m = %.2f K", T)

set label lbl_text at 60,1200

plot 'final.profile' u (\$2*${lx}):(\$4*${ke2T}) w p, 700 w l lw 2 lt 0, 1800 w l lw 2 lt 0, T w l lt 0 lw 2, g(x) w l lt 0 lw 2, f(x) w l lt 0 lw 2

set term post enha color 20
set out 'temp_profile.eps'
rep
EOF

# Now to analyze and draw the figures
gnuplot .plot.script
evince temp_profile.eps
#
rm -rf .plot.script fit.log
